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Abstract 

The standard approach for photoacoustic imaging with variable speed of 
sound is time reversal, which consists in solving a well-posed final-boundary 
value problem backwards in time. This paper investigates the iterative Landwe- 
ber regularization algorithm, where convergence is guaranteed by standard reg¬ 
ularization theory, notably also in cases of trapping sound speed or for short 
measurement times. We formulate and solve the direct and inverse problem on 
, what is common in standard photoacoustic imaging, but not for time-reversal 
algorithms, . We both the direct and adjoint photoacoustic operator an interior 
and an exterior equation . The prior is solved using a Galerkin scheme in space 
and finite difference discretization in time, while the latter boundary integral 
equation. We therefore use a BEM-FEM approach for numerical solution of the 
forward operators. We analyze this method, prove convergence, and provide 
numerical tests. Moreover, we compare the approach to time reversal. 

Keywords: Photoacoustic Imaging, Variable Sound Speed, Regularization 


Introduction 

Photoacoustic Imaging (PAI) is a novel imaging technique that allows for three dimen¬ 
sional imaging of small biological or medical specimens with high spatial resolution. 
It utilizes that an object expands and emits ultrasound waves when it is exposed to a 
short pulse of electromagnetic radiation (see e.g. |43l [40] ). The emitted ultrasound is 
assumed to be proportional to the electromagnetic absorption density which provides 
detailed anatomical and functional information. PAI aims for visualization of the 
absorption density by using measurements of the emitted wave outside of the object. 

Opposed to the standard PAI problem [24l |39l [40] we assume a spatially varying 
speed of sound. The underlying mathematical model of the wave propagation with 

*Laboratoire de Mathematiques LMIA, Universite de Haute Alsace, 4, rue des Freres Lumiere, 
68200 Mulhouse, France. (zakcLria.belhachini@uha.fr) 

^Computational Science Center, University of Vienna, Oskar-Morgenstern Platz 1, A-1090 Vienna, 
Austria, (thomas. glatzlOunivie. ac. at) 

^ Computational Science Center, University of Vienna, Oskar-Morgenstern Platz 1, A- 
1090 Vienna, Austria, and Johann Radon Institute for Computational and Applied Math¬ 
ematics (RICAM), Austrian Academy of Sciences, Altenbergerstrafie 69, A-4040 Linz, Aus¬ 
tria.(otmar .scherzerOunivie.ac.at) 


1 



spatially varying sound speed c(x) is the acoustic wave equation 


t) - Ay{x, t) = 0 in R” X (0, oo), 

c^{x) 

y( 2 ;, 0) =/(x) in R" , 
j/'(x, 0) = 0 in R" 

where / denotes the absorption density, which is proportional to the absorbed elec¬ 
tromagnetic energy. 

The inverse problem of photoacoustics with variable wave speed consists in de¬ 
termining the function / from measurement data m of y on a surface over time 
(0,T). That is, it consists in solving the equation 

y\T.=m. (2) 

With constant sound speed there exists a variety of analytical reconstruction formulas 
and numerical inversion techniques. We mention gSHHIHlllIlllS], see also the sur¬ 
veys [MllinillSlISHl- In the case of inhomogeneous sound speed, exact reconstruction 
formulas have been derived in [S]. A numerical method, providing approximations, 
is time reversal [MlIlTllIllEI], . Thereby the measurement data serve as Dirichlet 
boundary data, and the initial conditions at final observation time T are assumed to 
be identically zero. On top of the original algorithm of Fink m, the approach in |35j 
suggests a harmonic extension of the boundary data at time T as initialization data. 
The reconstruction obtained by time-reversal, lets say /t, approximates the true so¬ 
lution / as T increases to infinity. The exact / can be reconstructed if y{-,T) = 0, 
which happens, aside from trivial cases, only in odd dimensions and for homogeneous 
sound speed. In all other cases the results are error-prone. Moreover, this method 
produces approximations of / only under non-trapping conditions on c (see |21l I20j b 

Stefanov and Uhlmann [3^ showed that if diam(r2) = 2To (where diam(r2) denotes 
the diameter of fl with respect to the Riemannian metric c~^dx), then an observation 
time T > To is sufficient for a unique reconstruction of /. However, for stability of 
the inverse problem one needs longer measurement times and a non-trapping speed of 
sound condition. In fact, if the measurement time T is larger than a certain threshold, 
depending on the longest geodesic in the metric induced by c, the algorithm presented 
in |35] provides a theoretically exact reconstruction in terms of a Neumann series 
that contains multiple, subsequent time reversal and forward propagation of the data 
term. A computational realization of this approach has been investigated in |32j . This 
algorithm serves as a benchmark for our proposed algorithm. 

The presented approach for PAI inversion with variable sound speed relies on 
linear regularization theory m- Specifically, we obtain regularized convergence to 
the minimum norm solution even for short measurement times. Moreover, we ob¬ 
tain a new reconstruction method that does not require an artificial cut-off of the 
measurement data, nor harmonic extension of the data at the final observation time 
T. 

For the numerical computations, we decouple the wave equation into an interior 
part (solved by finite element methods), and an exterior part (with homogeneous 
sound speed) that is rewritten in terms of a boundary integral formulation. We then 
solve the coupled BEM-FEM system numerically. Note that by this approach we use 
exact, non-reflecting boundary conditions mm and therefore avoid the necessity of 
a perfectly matched layer to deal with the cut-off outside the domain of interest. The 
results are compared to conventional time reversal and the Neumann-series approach. 
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The paper is organized as follows: In Section we formulate the direct photoa¬ 
coustic operator L in a suitable choice of function spaces, and derive the adjoint. 
In Section we give a short overview about Landweber iteration and review some 
regularization results regarding convergence and convergence rates in view of PAI 
reconstruction. Moreover, we discuss the relation to the Neumann series approach in 
[35]. In Section we formulate the transmission problem for the wave equation used 
for numerical computation of both the direct and the adjoint problem. We state the 
boundary relations used for taking into account the exterior domain. In addition, we 
briefly describe the used discretization. Finally, Section provides a comparison of 
our reconstruction algorithm with the state-of-the-art reconstruction by time reversal 
and the enhanced time reversal Neumann series method from [35] . 

All along this paper we use the following notation and abbreviations: 

Notation 1 Let LI be a non-empty, open, bounded and connected domain in R" with 
-boundary dLl. The vector n{x), with x G dLl, denotes the outward pointing unit 
normal vector. We use the following sets 

n+ := R"\n, n- ■.= n and E := X (0,r) . 

We use the following Hilbert spaces: 

• Lf{Ll) = {/9 S L^(R"') : p = 0 in R"\n}, with inner product 

(Pi,P2)L2(n) = / ^Pi{x)p2{x)dx . 


• For Ll = LI or R"; 

— Let Hq{LI) be the closure of differentiable functions on R” with compact 
support in LI, associated with the inner product 

(pi,P2)jji(n) ■Vp2{x)dx . 

— H^{Ll) denotes the standard Sobolev space with inner product 

{pi,P2)H^n) = j^Pi{x)p 2 {x)-\-Vpi{x) ■ Vp 2 {x)dx . 

• L^ {dLl) denotes the standard Hilbert space of square integrahle functions on dLl 
with inner product 

(01,02)^2(90)=/ Mx)Mx)dS{x) . 

Jon 

L^(S) denotes the standard Hilbert space of square integrahle functions on E 
with inner product 

(0i,02)i2(s) = / / (fi{x,t)(j) 2 {x,f)dS{x)dt . 

Jo Jao. 
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• The induced norms are denoted by |MI 2 , 2 (q), IMI_l 2 (s); IMI/i-i(n) IHIffi(n)’ 
respectively. 

For n bounded, ouFIqITI), the norms IHIji/ii-q) and ||•||^l(Q) are equivalent (see 
1^ Theorem 6.28]): 

C’o ||p||//i(n) ^ ||p||^fi(n) — llplli/i(n) for all p ^ Hq{Fl). (3) 

• The trace operator 70 : i/^(]R") —^ L^(dfl) restricts functions defined on R” 
onto dTl, respectively. This operator is the decomposition of the standard trace 
operator 

7 ; H\n) L^{dFl) 

and the restriction operator 

R : ^ H^iFl), 

and thus as a composition of two bounded operators m Theorem 5.22] bounded. 
We abbreviate the norm with 


C,:=\hoR\\. (4) 

Notation 2 The absorption density f and the sound speed (? are supposed to satisfy: 

• c G C^{n), satisfies 0 < Cmin < c{x) < Cmax and c is non-constant in fl. 

• Without loss of generality we assume that c = 1 in . 

• The absorption density function f G H^in) is compactly supported in ft: 
supp{f) c rj. 

For the sake of simplicity of notation we omit space and time arguments of functions 
whenever this is convenient. 


1 Direct Problem of Wave-Propagation 

We analyze the wave operator L mapping the absorption density / onto the solution 
y of the wave equation 0 restricted to E. That is 

L : H^in) ^ , f^y\^. (5) 


In the following we show that L is bounded. Let us write 

E{t) ■■= [ \iy' f+ \\7y\^dx . ( 6 ) 

jR" C 

Computing the derivative of E with respect to t and taking into account Q gives 
E'{t) = 2 f ^ y"y' - Ayy' dx = 0 . 


Consequently 

E{t) = EiO) = , (7) 
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which implies that 


[ {y'fdx< [ \iyTdx< ||/||ii(n) 

jR" jR" C 


and 


II2^(^)IIr 1(R") — ll/llRi(n) 

for every t € (0, T). 

Lemma 1.1 Let y be the solution of Q, then 

II2/WIIri(R") < ll/llRi(n) ’ ^ ^ (0.^) • 

with 

C{T) 

where Cq is defined in (Ill- 

Proof. First, we note that for arbitrary t € (0,T), it follows from Q that: 

J {y{x,t) - y{x,Q)f dx = j y'{x,i)di^ dx 

<t f [ ljj'{x,i)Y didx 

jR" Jo 

^2 


/R" Jo 

= t I {y'{xfi)Y dx di 

Jo JR" 

+2 II ^||2 


< r 


1 r 1 (£ 2 ) 


< T 


-2 II i -||2 




Because (a — b)^ > — b"^ it follows from ([^ that 

[ {y{X:t)f dx <2 [ {y{x,t)-y{x,0)fdx + 2[ {y{x,0)f dx 
jR’* JR" JR" 

<27^^ll/llR„bn)+2|l/lli.(n) 

< 2max{l,T2} |j/|l^i(s^) . 

Because / G Hq{^) it follows that 

\W)\\l^(n-) < ^inax{l,T 2 } . 

This, together with (§ shows that for all t G (0,T): 

lly(0llRl(R") < V ll/llRl(n) ■ 


( 8 ) 


(9) 


( 10 ) 


□ 


In the following we prove boundedness of L: 
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Theorem 1.2 The operator L : H^iVl) —>■ T^(S) is bounded and 


\L\\ <C^CiT)Vf . 


(11) 


Proof. For given / let y be the solution of 0). From Q it follows that the solution 
y of (10 1 is in for every t > 0. Thus from ([^ and (10) it follows that 


r [ y\t)dadt<C^^ f \\y(t)\\l^i:R^)dt<C^,C{TfT\ 

Jo Jon Jo 




which gives the assertion. 


Remark 1.3 (Injectivity of L) In order to obtain injectivity of L, we need T to be 
sufficiently large. To specify this, we define 


To := max (dist(a;, 912)) , ( 12 ) 

where dist(x, 912) is the distance of x to the closest point x' € 912 with respect to 
the Riemannian metric c~^dx (see also From IU5\. Thm. 2] it follows that if 

T > Tq, than L[f] = 0 implies f = 0 in (0, T) x M”. 


□ 

In the following we characterize the adjoint of L : ilo(12) —?► -9^(S) on a dense 
subset of Lf{Yf). Because we know from elementary functional analysis that L* : 
T^(S) —Ro( 12 ) is bounded, 

Definition 1.4 Let i be the embedding operator from HQ{n) to Then i* : 

LflTl) —>■ ilg(12) is the operator which maps a function if G LfiLl) onto the solution 
of the equation 

—Am = if in LI , u = 0 on 912 . 

That is 

i* = -A-\ (13) 

where A is the Laplace-operator with homogeneous Dirichlet boundary conditions. 


In the following we derive the adjoint L* : of the operator L, which is required for 
the implementation of the Landweber iteration below. 


Theorem 1.5 Fork G C°°((0,T) x 912) the adjoint of the operator L, defined in ([^, 
is given by 

L*[h]=i* oL*jj[h] (14) 


where 


L*D[h] 



(15) 


and z := zifi) is the weak solution of 


1 


z” -Az = 0 in ]R"\912 x (0, T) 


[ z ]=0 


z{T) = z'{T) = 0 OT R", 
dz 


(16) 


dn 


= h on 912 X (0, T) . 
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Here 


[z :=z'\s-z Is 


and 


dz 

dz+ 


dz 

dn 

dn 

s 

dn 


where z+ := 2;|a+x(o,T) and z := 2:|nx(o,T)- 


Proof. For h G C°°{{0,T) x dil) the existence of a weak solution of (42) is proven 
in the Appendix. Taking v = y where y denotes the solution of ([^ it follows that 


hL[f]dS{x)dt = / hydS{x)dt 
i ds 

Jn 

-h 

= - / Vf A 


> 

1 

rz'(o)ii 




fdx 


,-i 


= / Vi* 


z'iO) 


z'{0) 


J 


z'{0) 


\7fdx 


Him 


V f dx 


(17) 


□ 


2 Landweber Iteration for Solving the Inverse Prob¬ 
lem of Photoacoustics 

The photoacoustic imaging problem rewrites as the solution of the operator equation: 

L[f] = m . (18) 

If the null-space of L is non-trivial, then iterative regularization algorithms, in general, 
when m is an element of the range of L reconstruct the minimum norm solution 


/t = L^m] , 


(19) 


where denotes the Moore-Penrose inverse of L (see 
Penrose inverses). 


for a survey on Moore- 


We propose to use the Landweber’s iteration for solving (18), because it can be 


compared with time reversal methods, which are the standard references in this field. 
More efficient regularization algorithms are at hand US], but these are less intuitive 
to be compared with time reversal. 

In the following we review properties of the Landweber iteration in an abstract 
setting (see [TS] |S]). We use the same notation for the abstract operator and the 
photoacoustic operator and measurement data m, , respectively, in order to have 
direct connection. 
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2.1 Abstract Landweber Regularization 

Everything that is formulated below is based on the following assumption: 

Assumption 2.1 Let L : Hi —> H 2 be an operator between Hilbert spaces Hi and 
H 2 satisfying uj ||E||^ < 1 for some oj > 0. Moreover, we assume that data of m 
is available (typically considered as noisy data), whieh satisfy 

||m — < S . (20) 

Then the Landweber iteration reads as follows: 


/o:=0 and fi = fi_i - ujL*[L[f^_i] - m^], fc=l,2 ,.... 


( 21 ) 


In case 5 = 0, that is, if = m, then we write fk instead of /^. 

Let T > 1 be some fixed constant. The Landweber iteration is only iterated for 
k = 1,2,... as long as 


h -L[fk]\\H, > ■ 


( 22 ) 


The index, where (22) is violated for the first time is denoted by k^{5,m°). 

The following theorem shows that the Landweber iteration converges to the best- 
approximating solution: 

Theorem 2.2 Let m G H{L) (note that the range of L equals the domain of Lf). 

• Let 5 = 0, then the Landweber iterates (fk) (cf. (HI); converge to the p, i.e.. 


\fk-p\ 




0. In addition, we have the series expansion: 


P =Y,{I-u:L*Ly[L* 

J=o 


• For 6 > 0 and satisfying ||m —< 5 let fcf 
Then 

p6 -Ghi P ■ 


ka.{6, m°) as in (22). 


Moreover, if m ^ ^{Lf), then ||/fe||^^ —)■ 00 as k —>■ 00. 

In the following we prove properties of the wave-operator L, such that we can apply 
the general regularization results. 


2.2 Convergence of the Landweber Iteration for the Photoa¬ 
coustic Problem 

In the following we apply Theorem |2.2[ for the photoacoustic imaging problem. 

Corollary 2.3 Let u! < = 0{\/Tp and L : iLg(fl) —tv^(E) as in (§. 

Moreover, assume that 

Then the Landweber iterates satisfy: 

• 7/5 = 0, then 

00 

P = Y,{I-LMy[L*[m]]. 

k^O 
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• For S > 0, the Landweher iteration is terminated at k^{6) := k^{6, m^) according 
to (221. Then 

for S^O. 


• If T > Tq, the reconstruction is unique, and 
solution. 


converges to the unique 


Proof. First, we note that from (11) it follows that lo ||L|| < 1. Then, the first two 
items follow directly from Theorems 2.2| 


From the injectivity of L for T > Tq (see Remark 1.3) it follows that p = /, 
which implies unique reconstruction. □ 


2.3 Comparison with Time Reversal 

We compare our approach with different variants of time reversal. We formally define 
the time reversal operator: 


L[h]=zi;0), 


(23) 


where z is a solution of 

— z" — Az = 0 in fl X (0, T), 

z(r) = z'(r) = 0 in 11, (24) 

z = h on dfl X (0, T), 

The fundamental differences between L and L* are that L is defined for functions with 
support in Q and that therefore L* requires a transmission condition in its definition. 

Stefanov and Uhlmann |35j modified the time reversal approach in the following 
sense: Rather than assuming (in most cases unjustified) the initial data z{T) = 0, 
they used the harmonic extension of the data term h{s,T), for s S 912, as initial 
datum at T. That is, for 

—At/) = 0 in 12 , with (/>(•) = m(-, T) on 912 

the modified time-reversal operator 


L[h\ = z(-,0) 


(25) 


is defined by the solution of equation 

-^z" — I\z = 0 in 12 X (0, T), 

z[T)=f>, z'(r)=0inl2, (2^) 

z = h on 912 X (0, T) . 

They were able to show that under non-trapping conditions and for sufficiently large 
measurement time T, there exists a compact operator K : Ffo(12) —>■ i?Q(12) satisfying 
\\K\\ < 1, and 


LL = Id-K, 


(27) 
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Therefore, the initial condition / can be expanded into the Neumann series 


f = [m] 

3=0 


(28) 


By induction, it is easy to see that the m-th iterate can be written as 

fk = fk-i - L[L[fk-i] - m ], (29) 

where 

fk=^K^ [m] . 

3=0 


At this point, we emphasize on the structural similarity between (29) and (21). 


Remark 2.4 We emphasize that for time-reversal there is no theory on stopping in 
case of error-prone data, such as we have available for the Landweber iteration. 


3 Numerical realization of L and L* 


We solve ([^ and (16) with the same numerical framework. By changing the variable 
t —>■ T — t in (fl6|), both equations can be rewritten as the transmission problem: 


-Av = 0 in ]R”\ail x (0, T) 

w(0) = Wo 5 w'(0) = 0 in]R"\9n, 

dv 


(30) 


dn 


= P, [^^] = 0 


on S . 


where for ([^ p = 0 and vq = f G Hq{Q) and for (16) we have p=hG T^(S) and 

Wo = 0 . 

Let vff = Wq in fl”, = 0 in 11+, then w+ = w|f2± satisfy, respectively: 


— Aw± =0 in 11+ X (0, T) 


1 + 

-li 


w±(0 )=w+, w±(0)=0 inll+, 

together with the transmission conditions 

’ dv' 


dn 


= p and [w] = 0 on S . 


(31) 


(32) 


Let G denote the fundamental solution of the standard wave equation with = 1 
in R". It is defined in R" x R, and its explicit expression 




47r|tc| ’ 


n = 2, 
n = 3, 


with H denoting the Heaviside step function, and S{-) being the 3D Dirac delta 
distribution (see, e.g., [Hj). 
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The (retarded) single- and double- layer potentials for {x,t) G E are defined by 


V[(p]{x,t) := [ [ G{y-x,t-T)ip{y,T)dSydT, 

Jo Jon 

IC[i(]{x,t) ■■= [ [ 

Jo Jon 


* '■ dG{y — x,t — T) 


(33) 


'>p{y,T)dSydT . 


Because we assume that c = 1 outside of $7 and because we assume that dH. is a G^ 
boundary, it follows that satisfies (see [T5]l: 


1 


v'^{x, t) = —V 


dv+ 


dn 


(x, t) — JC [u"*"] (x, t) , for all {x, t) € E . (34) 


The transmission conditions imply that 


V 


+ 


, dv+ 

= V and —— = 
on 


—-h p on E . 

On 


Therefore the equation for v = v on = 17 can be rewritten as follows: 


— Av = 0 
u(0) = Vo , v'(0) = 0 


Iv + V 


dv 

dn 


+ IC[v] = 0 


in (17) X (0,r), 
in 17, 
on E . 


(35) 


The numerical 
multiplying by 


solution is based on a weak formulation of (35). Integrating over 17, 
w G H^{Q) and integration by parts of the first line of (35) gives 


dt"^ 


-^v{t),w 


+ (Vu(7),Vu;)^2(o) - 

L2(n) 



L'2(dn) 


= 0 . 


Additionally we introduce the unknown function A := ^ defined on E. We are 
therefore searching for a solution (u, A) G i/^(17) x i7“^/^(917) that satisfies for almost 
all t G (0, T) the system 


^ /^u(7),u;\ + (Vu(i), Vu;)i2(fj) - (A,w)i2(an) 

u(0) = vo,v'{0) 

^vit) + V[X + p\{t) +^H(7) 


0 

0 

0 


for all w G H^{n ), 

in 17, 
on E . 

(36) 


For detailed analysis of the discretization of equations of that type ( [M| we refer to 
m, since we closely follow their approach. 

Next, we discuss the time discretization of the integrals V[A] and K\v\. We consider 
a uniform time discretization of the interval (0,T) into N steps of length At = T/N, 
and defining the discrete time levels = nAt. Following [12], by using Lubich’s 
convolution quadrature formula [29], we approximate V[A] and lC[v] at the time steps 
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tn, n = 0,..., N hy 


n „ 

V[ip]{x,tn) - y\)ip(y,)dSy , 

J=Oan 

n „ 

/C[V'](a;,i„) « ^ |a; - . 

_n 


(37) 


where the coefficients and w'A satisfy 

Wni^t, \x - y\) = ^ (\^- y\’ 

7—n V 


1^0 

L-1 


7 (pexp(i 727 r/L)) \ 

A, ; 


exp(—i?T, 27 r/L), 


(At, lx - y\) = ^ ^ ^|x - ?/|, 7 '(pexp^i^ 27 r/L)) \ gj^p('_j^27r/L), 

^ 1=0 ^ t / 


(38) 


where 


dv 

{r, s) = Ko{rs), K’^{r, s) = -s75ri(rs)— , 


and Ko{-),Ki{-) are the second kind modified Bessel functions of , respectively. The 
function 7 is given by 7 ( 0 ) = 3/2 — 2z + l/2z^ and is the associated characteristic 
quotient of the used backward differentiation formula method of order two. For the 
involved constants we choose L = 2N and /3 = where e is the machine precision 

(see uniin] for more details). 

In the first equation in (36), the second time derivative is approximated using the 
second order central difference expression 

^ - 2^” + ’ 

where u" := u(., f„). The first time derivative occurring in the initial condition is also 
discretized by central differences, namely 


^ 0 ^ / 1 -l^ 

—V = ^ (v — V 

dt At ^ ' 


0{At^) . 


From that we can restate the initial conditions as 
{v°,w) = {v°,w) , 

= (^v^,v^ — ^Af^ ((Vu°, Vic) (A°,w)) , for all w S . 

Using this, the explicit Euler discretization of ([^ is stated as 


(39) 


^ - 2u” +u”-i),u;) + (Vu", Vw) - (A”,u;) = 0, for all w S H\n ), 


,,7l+l 


(f) + = 0 on E , 


(40) 


for all 1 < n < TV. 

In space fl is triangulated, and we use piecewise quadratic basis functions, supple¬ 
mented by one cubic function for v. The functions A and are discretized by the 
use of piecewise linear functions. With this ansatz a mass-lumped integration scheme 
can be robustly implemented, which is not the case for purely piecewise quadratic 
functions. The details and numerical analysis to this scheme can be found in [^. 
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4 Numerical experiments and results 

Here, we present numerical results for the Landweber reconstructions and compare 


it with standard time-reversal reconstructions and the Neumann-series-approach (28) 
- always computed with the same discretization. We concentrate in particular on 
the cases where classical time reversal techniques have their theoretical and practical 
drawbacks, that is, so-called trapping speed geometries and short measurement times. 


In these cases, neither time reversal nor the Neumann-series (28) provide theoretical 


convergence. As we have shown in Corollary |2.3[ the Landweber reconstruction con¬ 
verges to the least squares solution, regardless of the chosen measurement time. In 
what follows, we also show the practical applicability of our method in such cases. 

To create the data sets, we use the discretization of the coupled method introduced 
in Section The forward computations are performed in a circle of radius 1. For 
both simulation and reconstructions, the chosen finite element space is that of con¬ 
tinuous, piecewise quadratic functions in space. For the spatial discretization of the 
boundary terms, we take piecewise linear basis functions. Note that for convenience 
and to optimize the computational effort, we assume dfl to be a circle of Radius R. 
For the time discretization, wo choose the step size At = h/(15cniax) for the simula¬ 
tion and At^ = /ir/(I4ci„ax)- Here Cmax is the maximum speed of sound. The time 
reversal reconstructions are obtained by solving the initial boundary value problem 


(24) with homogeneous initial values. In the images, when we used the harmonic ex¬ 


tension time-reversal (26), the heading is Neumann. The Landweber reconstruction 
is performed by the scheme described in Corollary |2.3| The choices of c in both the 


trapping and non-trapping case have been taken from 


4.1 Non-trapping sound speed 



Figure I: 


Let’s first consider the case of the non-trapping speed 


c{x) 


1 -I- 0.2 sin(27ra:i) -1-0.1 cos(27ra;2), x G Hfi-a , 
1 outside Bp> . 


Note that here and below, the speed of sound is smoothed out near the boundaries to 
satisfy our smoothness assumptions. The time reversal method already gives reliable 
reconstructions if the measurement time T is sufficiently large. For the test example, 
a measurement time of T = 4To (see Remark 1.3) is enough to provide a quantita¬ 
tively reasonable reconstruction by time reversal. This is illustrated in Figure We 
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Speed of Sound; Max Dist. to Boundary = 0.88061 


Speed of Sound; Max Dist. to Boundary = 0.90083 




Figure 2: 


Time Reversal 




Neumann Series Iterate 5 


Landweber Iterate 20. 



Figure 3: Data: Non-trapping speed, T = ATq. 


Time Reversal Neumann Series Iterate 5 Landweber Iterate 20. 



Figure 4: Data: Non-trapping speed, T = 1.2To. 


therefore are interested in the case where the measurement time is shorter, e.g. T as 
near as possible at Tq. In Figure]^ we compare the time reversal reconstruction and 
the Neumann series approach with our method, using T = 1.2ro- The Landweber 
reconstruction is stopped at a suitable stage of iteration. In practice, the improve¬ 
ment is very large in the first steps, while it needs a lot of iterations to satisfy the 
discrepancy principle from Theorem |2.2[ 

The main differences are clearly visible to the naked eye: The time reversal recon¬ 
struction fails to compute the central point in the image and produces artifacts. These 
artifacts can be avoided by the use of the harmonic extension as in (24). We here 
present iterate j = 5 of the Neumann series (28). Also the Landweber reconstruction 
avoids to amplify the artifacts in the image center. Moreover, it delivers the correct 
quantitative values of the initial pressure, whereas time reversal underestimates these 
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Landweber Iterate 1. Landweber Iterate 10. Landweber Iterate 40. 



Figure 5: Landweber reconstruction at different iterations. Data: Non-trapping 
speed, T = 1.2To. 



values. However, the smoothing step naturally included in every Landweber iteration 
seems to make this approach more stable than the time reversal Neumann series. In 
fact, at least with our method of numerical wave propagation, we have to stop the 
Neumann series after 5-10 iterates before numerical errors are amplihed too much. 


Remark 4.1 (Convergence rates in practice) The main difference between the 
convergence rates of the Neumann series and the Landweber iteration lies in the di¬ 
vision by after every backpropagation step, indicated by (15). 

In Figure^^we show reconstructions using the non-trapping speed and measurement 
time T = 1.2To. The measured error Hj/^ — L/fc||^ 2 ( 2 ) decreasing till k = 50. 

However, the major visible improvements seem to occur within the first 20 iterations. 
We therefore use the picture of iteration number 20 for our practical comparisons with 
the other methods. 


4.2 Trapping sound speed 

In the second example, we want to deal with the trapping speed 


c{x) 


1-I-0.5 sin(—STrxi) cos( 37 rx), x € Br-s , 
1 outside Br . 


In this case, there are geodesics present that do not propagate singularities to the mea¬ 
surement surface within finite time. The Landweber approach is the only one that 
gives a theoretical convergence result in this case. In practice, we see that conven¬ 
tional time reversal, at least for T = 1.2Tq, fails to give a detailed reconstruction. The 
Neumann series approach and the Landweber iteration behave similarly, again with 
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Figure 6: Data: Trapping speed, T = dTo- 




Neumann Series Iterate 5 



Landweber Iterate 20. 


Figure 7: Data: Trapping speed, T = 2Tq. 


Time Reversal Neumann Series Iterate 5 Landweber Iterate 50. 



Figure 8: Data: Trapping speed, T = 1.2To. 


the advantage of the Neumann series giving faster convergence, whereas the Landwe¬ 
ber gives a regularized solution that seems to be more robust against numerical errors 
and noise (figures and . 
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Appendices 

A Well-posedness of the transmission problem 

Notation 3 • The space 

V ■-{u & u\n & H\V)] (41) 

is a Hilbert space with respect to the inner product 

{(j)l,(j) 2 )y = / (I)i{x)4>2{x) dx + / V4>i-V4>2dx. 

Ja JR" 

The completeness follows directly from the properties of the Sobolev spaces 
andH\n). 

The dual space is denoted by V'. 


Definition A.l For h € C°°{(0,T) x dfX). A weak solution z of (16) satisfies: 


z e L^{0, T; V),z' € L2(0, T; g 1 , 2 ^ 0 , T; V'), 


where V is defined in (41), together with 


• z(T) = z'(T) = 0. Note that and thus z € H^(0,T; L^(]R^)) and z' € 
H^{0,T;V'), such that the traces make sense. In fact z(T) G L^(R") and 
z'{T) G V. 


-t:z''v dxdt + 


/R" 


/R" 


Vz(t) • Wv{t) dx dt 


(42) 


/an 


h{t)v{t) dS{x) dt , for all V G L^{0,T-,V). 


Note, that the definition of a weak solution is not used consistently in the lit¬ 
erature. Evans m, for instance, uses the test functions v time independently, 
and formulates a weak form of (|16[) for almost all t G (0, T). 


Theorem A.2 For h G C°°{{0,T) x dU) there exists a weak solution z of (16) - that 
IS of (@. 

Proof. The proof is similar to [THl Section 7.2.2]. 

1. First we construct a Galerkin-approximation: Let be an orthonormal 

basis of L^(]R"), which is an orthogonal basis for For fixed m let 


Zm{t) = '^d'f,{t)wk , tG{0,T) 




where 


d^{T) = 0, dt'{T) = 0 meaning that Zm{T) = z'^{T) = 0 . (43) 
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and 


+ {zTmWk) hUb.^) — (^i'f^fe)L2(s) , 

l2(R") ' (44) 

for all 0 < t < T, for all fc = 1,..., to . 

Analogously to uni Theorem 1, Section 7.2.2] the Galerkin approximation can 
be shown. 



2. The following estimates are different to uni Theorem 1, Section 7.2.2] and thus 
included - it is essential to consider an additional time integration. As in m we 
use it)wk as a test function in the Galerkin approximation. 

However, we do not apply it pointwise for every r S (0,T), but in integrated 


form. Then from (44) it follows that 

f f ^mii) +Vzm- Vz'^dxdi = - f f h{i)z'^{i)didS{x) 

for all T G (0, T) . 

Partial integration of the right hand side and evaluation of the integral terms 
on the left hand gives 

2 \ 


d / 

Z'raii) 

Jr di \ 

C2 


M 


L2(R") 


R1(R") 


dt 


(45) 


=2 / / h'{i)Zm{i) didS{x) — / h{T)zm{T)dS{x) 

\^JdQ JT JdQ ^ 

Since Zm{T) = z'^{T) = 0, the left hand side equals 




lkm('^)llRi(R")) • 

Estimating the right hand side by Cauchy-Schwarz-inequality and using mean 
inequality, we get for an arbitrary D > 0: 

2 


z'mi-r) 


L2(R") 


+ lkm('r)l|^l(pT.) 


< 


^ ^ di+ j 

11 (i ) 11 ^2 (an) + 11 ('^) 11 L2 (SO) 


+ D^ 


Let 


pT 

C{h,T):= J \\h'{i)\\l^^g^^di+\\h{T)\\l2^g^^ , 


C{h,t) := J C{h,T)T, 

and using the trace theorem Q it follows that 


Zm{'r) 


L2(R") 


+ \\Zm{^ 


Ir1(R") 


<^C{h,T)+D^C^^ 


di + \\Zm{T)\\ 1 ■ 
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By integrating of r over [t,T] we get the following estimate for all t S (0,T): 


r 

z'mi-r) 

L 

C2 


L2(R" 


dr 


<—C{h,t) + D^C^^ 


T 



|pm (0 ll//i(n) J 

In addition, because Zm(T) = 0, we know that 

Zm{T)dx] = / Zmir) - ZmiT)dx 



Moreover, we know from the inequality from that 


Zmir) - / Zm{T)dx 
JQ 


L2(n) 


^ C'g ll■^m(7■)||Rl(Rn) J 


which implies that for all t G (0,T): 

INm('^)llL2(n) — Zm(T)dx'j + 2C(3 || (t) || , 


and consequently 


l^m('r)|l^l(-Q) < Zm(T)dx'j + (2 Cq + 1) || (t) || 


Thus it follows from (46) that 
2 


r 

z'miT) 

L 

C2 


L2(R") 


IZm(T)llffl(j^^)dT 


T 



(46) 


(47) 




+ D (7^(2(7(3 + 1) / I / |pm(i)||+ ||zm(T)||p 3 i(p{Ti) I dr 


Zm{T)dx ) dr . 

0 \Jn 


(48) 
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Taking t = 0 then gives 


r 

z'miT) 

Jo 

C2 


L2(E") 


+ ||■^m(7 


IhI(h-) 


dr 




/•T 

+ D^C^{2Cg + 1){T+1) \\zm{i)\\l, 

Jo ° 

+ J (^J Zjn{T)dx^ dr . 


(R") 


dt 


Thus it follows from (47) that 




z'mi'r) 


L2(R") 


+ INm(T)||^l(R„)dr 




+ D^C^^{2Cg + 1){T+1) f \\zmii)\\l, 

Jo “ 

+ J (^J ZmiT)dx^ dr . 


(R") 


dt 


Choosing D such that 


max {Z7"C2(2 Cg + l)(r + l),2D^C^{T + 1)} < ^ . 


provides that 


J Zm(T)dx^ dT+ 


1 

4r 
1 
4 


r 

z'mi'r) 

Jo 

C2 


L2(R") 


+ INm(T)|l^l(R„)dr 


sup 


(49) 


(50) 


Therefore there exists a constant C such that 

max ||| 2 rnllL 2 (Q_y.j^n) , || 2 m|li 2 (o_T;^fi(p{np , || | Q || ^2 .//l (q)) | < C , (51) 

where the last estimate is again due to inequality , already used above. 
Because {wk} is a basis it follows from (441 that for all v G L^(0,T;V): 

J! 


L2(0,T;R^ 


< sup (Zm, ’^)r1(0,T;R") sup (/l, v) i2(2) — ^) : 


where C{h,C) is constant depending on h and C, but is not dependent on m. 
This shows that G L'^{0,T; V). 
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3. (50) guarantees that {zm} is uniformly bounded in L^(0, T; iJg (H")), and 
{• 2 m|n} is uniformly bounded in L^(0, T; i/^(n)). Moreover, {z'^j is uniformly 
bounded in L^(0, T; L^(R")), and {z'^} is uniformly bounded in L'^{Q,T-,V'). 
Thus it has a weakly convergent subsequence, which we again denote by {zm} 
which is converging weakly in the following sense: 


{zm} z with respect to L^(0,T;_ffg(lR”)), 
{zm\Q} z\q, with respect to L^(0, T; , 

{z'm.} '0 respect to L^(0,T; L^(1R")), 

{z'm} P with respect to L^{0,T\V') . 


We note that the trace operator 7 : H^{n) —)• L'^{dfl) is bounded. Therefore, 
every test function v G LF'(Q,T\V) satisfies u|an G L'^{0,T; L^{dfl)). Conse¬ 
quently, 


/o Jon 


ZmV dS{x), dt 


lo Jon 


zv dS{x), dt . 


Now, let V G (7“((0, T) x R"), then v', v" G T) x R") as well, such that 


/ ipvdxdt = lim / / z'^vdxdt 

Jq, Jo Jq 


= — lim 


Zmv'dxdt 


Jo Jn 


= lim 


= lim 


zv'dxdt 


lo Jn 

pT p 


z'vdxdt, 


lo Jn 


which implies that tjj = z'. Analogously one can show that 
L^(0,T;F'). Therefore, z it is a solution of (42). 


□ 


Theorem A.3 The solution of (16) is unique in L^{0,T;V). 


Proof. Let us assume that there exist two solutions zi,Z 2 of (16). Then 


4 (zi - Z 2 )" - A(zi - Z 2 ) = 0 in R”\9r! x (0, T ), 

(zi - Z2)(r) = (zi - Z2)'(r) = 0 in R", 

'd{zi - Z 2 )' 


(52) 


[{Zi - Z2)] = 0, 


dn 


= 0ondnx (0,r) . 


But this solution has only a trivial solution, which one sees from ([^ and by noting 
the boundary conditions at 00 of an iJg(R") function. □ 


24 









